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Abstract — We present a study about the possibility to detect 
neutrino induced extensive air showers at the Pierre Auger 
Observatory. The Monte Carlo simulations performed take into 
account the details of the neutrino propagation inside the Earth, 
the air as well as the surrounding mountains which are modelled 
by a digital elevation map. Details on the sensitivity with respect 
to the incoming direction as well as the aperture and the total 
observable event rates are calculated on the basis of various 
assumptions of the incoming neutrino flux. 

I. Introduction 

The detection of very high energy cosmic neutrinos, above 
1 EeV, is important as it may allow to identify the most 
powerful sources in the Universe. 

In the first place, neutrinos, due to their small cross-section 
can travel cosmological distances without interactions. Hence 
they might carry astrophysical information about their sources 
which are commonly believed to be optically thick [1]. 

Second, due to their connection to the emission of cosmic 
nuclei and gamma rays, they might help to solve the problem 
of the origin of ultra high energy cosmic rays (UHECRs). 
In fact although the existence of UHECRs is experimentally 
proven, their composition and origin are still unknown. Many 
models have been proposed to explain the origin of UHE- 
CRs. Some of these, which involve mechanisms of particle 
acceleration (bottom-up), claim that they might be produced 
by Active Galactic Nuclei (AGN) and Gamma Ray Bursts 
(GRB) [2]. In such scenarios neutrinos could be produced 
with an upper flux limit given by the Waxman-Bahcall (WB) 
bound [1]. Other models claim that UHECRs might come 
from the decay of super-massive objects (top-down): These 
objects are expected to be produced by radiation, interaction 
or collapse of topological defects such as monopoles, cosmic 
strings, etc. [3]. The topological defect models predict a larger 
flux of photons and neutrinos arriving at Earth than the bottom- 
up models. Finally, high energy neutrinos could be also 
produced through pion decay when protons interact with the 
cosmic microwave background. These neutrinos are so-called 
cosmogenic neutrinos [4]. 

Nevertheless, due to their vacuum oscillations [5], a flux of 



high energy cosmic neutrinos is expected to be almost equally 
distributed among the three neutrino flavours. Therefore the 
study of oscillation effects on high energy neutrino fluxes can 
be used to study the neutrino mixing and distinguish amongst 
different mass schemes. 

According to the models described above, a low incoming 
flux of neutrinos is expected. In addition, neutrinos have a 
small interaction probability. Therefore, in order to get a 
detectable rate, very large neutrino detectors are needed such 
as the Pierre Auger Observatory [6]. In particular, for Auger 
Observatory Earth-skimming v T showers are the best suited 
candidates to produce detectable showers. Studies of such 
possibilities were recently presented in Refs. [7]— [ 11]. 

In this paper to exemplify the potential and the features 
of the developed simulation tool we study neutrino induced 
extensive air showers (EAS) focusing on the fluorescence 
detector (FD) set up at the Pierre Auger Observatory. Taking 
into account the details of the neutrino propagation inside the 
Earth and in the Atmosphere as well as the topography of 
the Pierre Auger Observatory site, the directional dependence 
of the neutrino rate for up-going and down-going neutrino 
showers (probability maps) are evaluated accounting for the 
10% duty cycle for the FD. The WB bound is used as an initial 
neutrino flux. We use the ellipsoidal model of the Earth based 
on the so-called datum WGS84 [12] model. In addition, we 
use an elevation map of the mountains surrounding the Auger 
site. We find that the shape of the surrounding area influences 
significantly the calculated neutrino rates. To investigate the 
response of the FD for up-going showers, we generate the 
longitudinal profiles of shower development using the EAS 
MC generator AIRES [13]. The light propagation and the 
hardware detector trigger are simulated by means of the Auger 
software framework called Offline [14]. Finally the aperture, 
acceptance and neutrino rate are calculated for up-going v T 
showers. 

II. Method 

The software tool used in the present paper is based on the 
code ANIS (All Neutrino Interaction Simulation) [15]. ANIS is 



Fig. 1. Topography of the Auger site according to CGIAR-CSI data. The 
centre of the map corresponds to the centre of the Auger array (latitude 

<t>Anger.center = 35.25° S, longitude ^Auger.center = 69.25° W). The 

Auger position is marked by a circle. 

a code originally developed by A. Gazizov and M. P. Kowalski 
for the complete simulation chain of neutrino propagation and 
interaction in the context of the AMANDA experiment [16]. 
It allows to generate z/-events of all flavours, to propagate 
them through the Earth and finally to simulate ^—interactions 
within a specified volume inside the Earth. All relevant Stan- 
dard Model processes (charged current (CC), neutral current 
(NC) ^^-interactions and resonant v e — e~ scattering) are 
implemented. In addition, neutrino regeneration at all orders is 
included. The density profile of the Earth is chosen according 
to the Preliminary Earth Model [17]. Deep inelastic v — N- 
cross-sections are calculated from Quantum Cromo Dynamics 
(pQCD) with structure functions according to CTEQ5 [18] 
and with a logarithmic extension into the small-x region. The 
tau decay is simulated using the program TAUOLA [19]. In 
order to simulate neutrino showers suitable to be detected by 
the Pierre Auger Observatory, some important changes and 
extensions of the code were needed. For instance, in ANIS 
neutrino propagation and tau decay are simulated only inside 
the Earth. Instead, we need to simulate also the propagation 
through the atmosphere. Moreover, it is necessary to take 
into account that the Auger Observatory is positioned on the 
surface of the Earth. 

First, the topography of the Auger site was implemented. 
The description of the relief of the Andes mountains was 
made according to a digital elevation map (DEM). These data 
are available from the Consortium for Spatial Information 
(CGIAR-CSI) [20]. The implemented map of the area around 
the Auger site obtained is shown in Fig. \l\ 

Second, a redefinition of the detection volume was done. 
The original ANIS version uses the concept of the detection 
volume [21]. This concept is useful for calculating event rates 
for a given neutrino flux §{E U ,9 V ). The detection volume 
corresponds to the so called active volume in which potentially 
detectable neutrino interactions are simulated. In the original 
ANIS this volume is defined as a cylinder with z-axis parallel 
to the neutrino direction [15]. In our adapted version we kept 



Fig. 2. Sketch of the geometry relevant for neutrino simulations, see text 
for more details. 

the general idea of active volume but with some modifications, 
as shown in Fig. [2] As it can be seen from this figure, the 
active volume for a given incoming neutrino with energy E v is 
defined by a particular plane A gen and distance AL. The plane 
A gen is the cross-sectional area of the detector volume and it 
was used as reference plane for the generation of incoming 
neutrinos. The area depends on the zenith angle 9 of the 
incoming neutrino. The detector was modelled as a cylinder 
with radius R and height H, The distance AL is the multiple, 
n, of the average lepton range (Ri ep {Ei ep )) [22]. Inside the 
active volume lepton tracks are also shown. Neutrinos with 
energy E u produce leptons with individual ranges, depending 
on the fraction of energy transfered. More precisely, the 
incoming neutrino is forced to interact in the active volume 
according to its interaction probability 

P(E V , E lep , 9)~N A x a(E u ) x p(Z) x AL, (1) 

where cr{E u ) is the total neutrino cross section, p(Z) 
the local medium density and Na the Avogadro constant. 
P(E„, Ei ep ,9) is the probability that a neutrino with energy 
E v crossing the distance AL would produce a lepton with an 
energy Ei ep . In this way the production vertex of the lepton is 
created. Then, the lepton propagates through the matter, losses 
some energy and decays (if it's unstable) in the vicinity of the 
detector or inside the detector volume. In order to calculate 
physical quantities, one has to weight the events. A first weight 
is the interaction probability defined by Eq. 0. A second 
weight comes from the normalisation of the injected neutrino 
flux $(E U , 6 V ) 

K = N'lxATx / <S>(E v )dE v x / A gen (9)xdn, 
JE min Je min 

(2) 

where dVl — 2ir sin 9d9 is the space angle, AT the observation 
time and N gen is the number of generated events from 
surface A gen . Here, we assume the isotropic neutrino flux 
Q(E V ,9 V ) = $(£?„). The expected event rate of the neutrinos 
in the detector volume can be calculated from 

Nratio = K X ^ Pi, (3) 
i=l 



where N acc is the number of events triggering the detector 
and passing all quality cuts of the cascade analysis. The 
normalisation factor F™ is chosen such that the rate of the 
neutrinos results in the total number of events per year. 

Third, new tau lepton decay channels were implemented. 
In the original ANIS code the tau lepton decay is simulated 
using the code of TAUOLA [19]. The tau lepton decay is 
simulated by randomly choosing one event from a data base 
of pre-simulated events. 




Fig. 3. 

energy. 



Tau lepton range R T and its decay length d T as a function of the 



— r T = -a - (3(E T ) x E T , 



(4) 



The lepton range, R T , is given by the integral of the inverse 
of the tau lepton loss rate over the tau lepton energy 

1 f 1 



Finally, the tau lepton decay routine was partly rewritten 
to take into account the processes of tauons escaping from 
the Earth's surface and decaying in the atmosphere. In ANIS 
the tau is propagated in small energy steps until the age of 
the tau lepton exceeds the tau lepton lifetime. This procedure 
works quite well if the production and decay vertex are in 
the same medium (either in air or rock), but if the tau lepton 
crosses the Earth's surfaces during its propagation inside the 
active volume, different amount of energy loss in the Earth's 
crust and air have to be taken into account. In particular, when 
a tau lepton is generated in the Earth, it loses energy due to 
ionisation and radiation processes. These energy losses per unit 
length of crossed matter (in g/cm 2 ) are usually approximated 
by a linear equation (continuous energy loss approach), which 
reads as 

dE T 
~dX 

where the factor a is due to the ionisation losses and (3 is due 
to radiation losses. The factor a is negligible for ultra-high 
energies. The factor (3 parameterises the tau lepton energy loss 
through bremsstrahlung, pair production, and photonuclear 
interaction. 

Exemplary we have used two parameterisations: j3 is linear 
dependent on energy 

f3 A {E T ) = 0.71xl0" 6 cm 2 g _1 +0.35xl0" 18 £; cm 2 g" 1 GeV" 1 

(5) 

used by Aramo et al. in Ref. [8], and 

Pb(E t ) = (1.508 + 6.3(£: T /10 9 ) - 2 ) x 10" 7 cm 2 g- 1 (6) 
from [15], [22], [23]. 



R T = 



T dE. 



(7) 



p{z) J dE/dX 

Thus differences in (3 will lead to different tau lepton ranges 
(length of the lepton track) and consequently influence the 
expected neutrino event rate and aperture calculations. 

During a single step the distance passed by the tau lepton 
is evaluated according to the following formula 

1 



AR T = 



-ln{El/Ei), 



(8) 



p(z)(3(E() 

where E\ is the initial energy of the tau lepton and E' T is 
the energy of the tau lepton after the distance AR T . During 
the energy step (3 is constant. The distance AR T is evaluated 
with the altitude dependent density p(z) = p rock (z) or p(z) = 
p air (z). p rock is the density of the Earth according to the 
Preliminary Earth Model [17] and p air is the air density 
calculated according to the US standard atmosphere (Linsley 
parameterisation) [24]. 

In Fig. |3] the tau lepton range in standard rock, p — 
2.6 g/cm 3 is shown for the parameter set (3b based on 
calculations reported in [23]. One can see that the tau lepton 
range R T is of the order of about 10 km at lEeV. This 
value is about 5 times smaller than the tau lepton decay 
length, d T , at the same energy. The decay length is given by 
d r = cT T (E T /m T ) ~ 49km x (£ T /10 18 eV) for tau lepton, 
where ct t = 87.11 is the tau lepton lifetime in [im and 
m T = 1777.03 MeV is the tau lepton mass [25]. Thus, a tau 
lepton with an energy of about 1 EeV produced close to the 
Earth's surface can escape from the Earth before decaying. 
One can note also, that the tau lepton decay length is about 
490 km for a tau lepton energy of about 10 EeV. Therefore an 
emerging tau lepton can escape from the atmosphere, which 
is assumed to have an height of 100 km. 

Applying the neutrino generation code, we obtain the tau 
lepton decay vertex position, the energy and momentum of the 
decay products for simulated neutrino showers with a given 
energy. AIRES was used to generate the longitudinal profiles 
of charged particles and the energy deposit based on the ANIS 
output. A special mode was used to inject simultaneously 
several particles at a given interaction point. The emission of 
fluorescence light by the shower together with its propagation 
towards the detector and the response of detector itself, in- 
cluding electronics and trigger, were simulated by the Offline 
software. Finally on the basis of the algorithms implemented 
in the Offline software the first two trigger levels called First 
Level Trigger (FLT) and Second Level Trigger (SET) were 
simulated [6]. After the simulated events had passed the FLT 
threshold trigger (at least one pixel), a search for pattern 
consisting of 5 pixels was performed according to the SET 
algorithm. In this way the trigger efficiency of the FD for a 
given neutrino energy can be defined as 

E(i5„) = x 7 (9) 

-L* Aires 



where N^ires is the the number of AIRES tau showers 
simulated for a given neutrino energy, N$lt the number of 
showers passing the SLT condition and 7 the duty cycle of 
fluorescence detector. 
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Fig. 4. Event rate as a function of azimuth and zenith for the detector 
volume in the case of quasi-horizontal u T showers, assuming E\ h = 1 EeV, 
a WB flux, AL = 15 (R T {E T )), B , 10% duty cycle of FD detector. 

III. Results 

Simulations are done for neutrino showers ranging from 1- 
100 EeV. Usually 500.000 events are generated on the surface 
A gen for different azimuth and zenith angles. A cuboid with 
an area of 50 x 60 km 2 and height of 10 km positioned at 
1430m a.s.l. was used as detector volume (= Vfd)- It agrees 
quite well with the detection volume seen by FD. 

In Fig. |4]the map of the expected event rate in yr -1 as a 
function of the incoming direction (9, 4>) is shown. Significant 
directional differences in the number of expected events are 
seen. The largest rate is for almost horizontal neutrinos with 
the zenith 9 in the range between 89°-93° and azimuth (f> in the 
range between 100° and 270° (North- West-South). However, 
some peaks in the rate distribution exist for azimuth of about 
350° (East). This behaviour is due to the fact that the Auger 
site is surrounded by the Andes on the West and smaller but 
closer mountains in the East (see the map shown in Fig. [Q. 
On the other words the correlations between the calculated 
distribution of the event rate and the topography of the Auger 
site is well seen (mountain effect). The mountain effect is 
much more pronounced for v T showers very close to the 
horizon and less obvious for larger zenith angle ranges. A 
quantitative analysis of this effect for horizontal showers is 
presented in Tab. |I] and for up-going showers in Tab. |H] 

In Tab. U the dependence of the expected number of events 
from the direction of incoming neutrinos is well seen. For 
example, in the case of quasi-horizontal (QH) showers the 
expected ratio from the West is about 3 times larger than the 
ratio from the East (for downward-going shower (DW) the 
ratio from the West is about 5 times larger than from the East). 
Note also that the rate in the case of QH showers is about one 
order larger than the rate for downward-going showers. 



TABLE I 

Expected event rate in yr -1 for quasi-horizontal showers. 
The rate has been evalueated for different incoming u t 
directions and at e* h = 1 eev assuming the fd detection 
volume V PD . The wb neutrino flux, j3 B and 10% duty cycle of 

FD DETECTOR HAVE BEEN ASSUMED. 





Azimuth 

(deg) 


ly nw 
(yr" 1 ) 


(yr" 1 ) 


Theta 




85° -90° 


85° -95° 


East 
North 
West 
South 


-45-45 
45-135 
135-225 
225-315 


0.0005 
0.0007 
0.0027 
0.0007 


0.009 
0.012 
0.023 
0.009 


Total 




0.0046 


0.052 



TABLE II 

Expected event rate in yr -1 for up-going neutrino showers 
for different ranges of zenith angle for e t h = 1 eev, 0b, fd 
detector volume v f d, wb flux and 10% duty cycle of fd 
detector. n f ^ is the rate calculated according to the 
spherical model of the earth with the detector positioned at 
10m above the surface. 



class A+B+C class A Sphere 

6 N FD N acc N FD N acc 3^ 

(deg) (yr- 1 ) (yr' 1 ) (yr' 1 ) (%) 

W92 0.0219±1% 34205 0.0108 1209 0.0105 48" 

92-94 0.0168±2% 9895 0.0133 1248 0.0140 24 

94-96 0.0068±15% 199 0.0060 165 0.0059 11 

96-98 0.0008±43% 29 0.0008 27 0.0010 6 



In Tab. |ll] the expected event rates for different ranges of 
zenith angle are listed. The rate is dominated by events close to 
the horizon and gets only small contributions at larger zenith 
angles. This is due to the fact that in the case of up-going v T 
the regeneration effect plays a important role. A high-energy 
v T interacts in the Earth producing a tau lepton which in turn 
decay into a v T with lower energy due to its short lifetime. 
The regeneration chain v T — > r — > v T — > ... continues until the 
tau lepton reaches the detector (in our case the active volume). 
This effect leads to a significant enhancement of the tau lepton 
flux up to about 40% more than the initial cosmic flux of tau 
neutrinos of energies between 10 6 — 10 8 GeV [26], [27]. 

Additionally, the rate for the different classes of events 
defined in Fig. |2] is presented. Class A has the production 
vertex in the Earth (below horizon) and decay vertices inside 
Vfd- The class B and C have production vertices above 
the horizon and decay vertices inside Vfd (the contribution 
of class C is about one order smaller than class B). In 
principle the difference between these classes are a measure 
of the influence of the mountain effect reflected in a zenith 
dependent event rate [28]. As it can be seen from this Table 
(the last column), the contribution of mountains to the total 
v T event rate is significant. For example, for nearly horizontal 
u T showers within zenith angle less than 2 degrees below the 
horizon, the contribution is about 50% while for larger zenith 



angle ranges it is on average less than 24%. We come to the 
same conclusion if we estimate the mountain effect taking into 
account the calculation with the simple spherical model of the 
Earth because the rate calculated in this case agrees very well 
with the rate calculated for class A. 

In Fig. |5J\ the estimation of the aperture for up-going u T 
showers is shown. We define the aperture by 

A{E V ) = N- e l xJ2H Pi,j( E u,E T ,6) x Aj(9) x Ail, 

i=l j=l 

(10) 

where N T is the number of tau leptons with energy E T . E u is 
the initial neutrino energy, N$ c £ is the number of tau leptons 
coming from a given direction inside Vfd and passing our 
cuts. The aperture was calculated according to the parameter- 
isation of 0a and /3b. Different parameterisations lead only 
to small differences in the calculated aperture. Additionally, 
in Fig.|5j\ the aperture calculated by [11] is shown. One can 
see a fairly good agreement within 20% between the results 
presented in this paper and aperture calculated by [11] for 
energies between 1 EeV and lOEeV. For larger energies, above 
lOOEeV, our aperture decreases with increasing of neutrino 
energies. This is due to the fact that for these energies the tau 
lepton decay length in air is extremely large (above 5000 km), 
then the emerging tau lepton from the Earth will escape from 
the detector, or even from the atmosphere. 
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Fig. 5. (A) The effective aperture for the Pierre Auger Observatory, A{E Vt ), 
for E T > 10 18 eV, 90° < 6 < 95° ,0° < <p < 360°, V FD and 
N gen = 200000; (B) The effective acceptance A C (E^) for the Pierre Auger 
Observatory. The 10% duty cycle of FD was included. 
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Fig. 6. Tau neutrino and anti-neutrino fluxes from various models, taken 
from Ref. [8]. 

TABLE III 

Expected event rate in yr _1 for the different neutrino fluxes 

based on our aperture (n ape.r) and acceptance calcultion 
(Nacc)- Results for u t events with zenith angle in the range 
90° - 95°, assuming^ = 1 EeV, a WB flux, AL = 15(R T (E T )), 

f3 B , 10% DUTY CYCLE OF FD DETECTOR. 





WB 


GZK-L 


GZK-H 


TD 


NH 


NAper 
Nacc 


0.044 
0.004 


0.456 
0.042 


0.699 
0.072 


0.546 
0.053 


1.478 
0.137 



In Fig. 133 the FD acceptance for Pierre Auger Observatory 
is shown. The acceptance is calculated from 

A CC {E„) = A{E„) x E(£„)| Br >xEeV. (11) 

assuming that a trrigered event is seen at least by one flu- 
orescence eye of Pierre Auger Observatory. The acceptance 
increases from 0.001 km 2 sr yr to about 0.003 km 2 sryr be- 
tween 30 EeV and 1000 EeV. The acceptance does not saturate 
but drops at larger energies since tau leptons at such energies 
escape unseen from the detector volume Vfd lowering the 
acceptance. The calculated acceptance differs about two orders 
of magnitude from the calculated aperture. 

Finally, in Tab. |lll] the rate for different injected neu- 
trino fluxes shown in Fig. |6] based on our acceptance 
and apperture calculation, are listed. The WB rate is ob- 
tained for Waxman-Bahcall bound. In this paper we use 
the conservative estimate of this bound ^(E^ T+ p T ) = 1 x 
10 -8 i?~ 2 (GeV s _1 cm -2 sr). Other rates are reported for 
the same neutrino fluxes considered in section 2 of Ref. [8]. 
The two GZK fluxes refer to the two possible scenarious 
for cosmogenic neutrinos, which are those produced from an 
initial flux of UHE protons. Instead the NH (New Hadrons) 
and TD (Topological Defects) cases are two examples of exotic 
models. As one can see from Tab. Mil the number of expected 
events equals about 1 neutrino event per 2 years in the case of 
our aperture calculation. Taking into account the FD detector 
properties rate is about one order smaller, but neutrino event 
still can be observable by Auger detector during a few decade 
of operational. 

Note that event rate becomes larger for lower threshold 



energy. In fact FD detector can record also the shower with 
enegies above 10 17 eV. 

IV. Summary 

A tool to simulate neutrino propagation and interaction in 
the Earth and in the atmosphere taking into account the local 
topographic conditions, was presented. A detailed investigation 
of neutrino induced showers was performed. The focus was 
put on the neutrino sensitivity for the FD detector of the Pierre 
Auger Observatory and therefore we have used the digital 
elevation map for this site. The probability map (event rate 
distribution) and the aperture were calculated and the impact 
of the mountains surrounding the Auger site was quantified. 
General features are consistent with previous results, but as 
a consequence of the surrounding mountains a pronounced 
variation of the event rate as a function of the azimuth is 
seen. 
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